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Abstract: 
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performance on many applications indicates their potential for real-time simulations. 
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1 Introduction 



Intersection problems are among the fundamental topics in computational geometry 
and geometric modeling. They are also considered important in computer animation, 
physical based modeling, robotics, and computer simulated or virtual environments. In 
most of these applications, the objects are undergoing motion governed by dynamic 
constraints or external forces. It is important to determine contact or interference 
between any pair of objects in a dynamic environment and generate a collision response 
to it. This problem has received a lot of attention over the last few years in all these 
areas. 

The problems of contact determination and interference detection have been ex- 
tensively studied in different areas. The literature in computational geometry consists 
of a number of theoretically efficient algorithms for polyhedral objects in static envi- 
ronments [1, 2, 8, 9, 10, 11, 12, 16, 34]. There are a number of algorithms with good 
asymptotic bounds, however their practical utility is not clear since many of them are 
not implemented in a realistic environment. 

Many algorithms are known for intersection between curved objects represented 
as algebraic surfaces or piecewise spline models in geometric modeling [13, 24, 31, 
23, 35]. However, they are targeted towards robust computations of the intersection 
curve between such models for CSG operations and boundary computations. They 
are usually specialized to static environments and are relatively slow for the kind of 
performance desired in computer simulated environments 

Similarly, the problem of collision detection has been addressed in robotics liter- 
ature as well. However, the main purpose of applications has been on planning a 
collision-free trajectory between obstacles [7, 14, 26, 29, 40]. This is different from 
applications in physical based modeling and virtual environments, where the motion 
is subject to dynamic constraints or external forces and cannot typically be expressed 
as a closed form function of time. . 

The simplest algorithms for interference detection in dynamic environments are 
based upon using bounding volumes and spatial decomposition techniques. Typical 
examples of bounding volumes include cuboids, spheres, octrees etc. and they are 
chosen due to the simplicity of finding interference between two such volumes. These 
bounding volumes work very well in performing “rejection tests”, when two objects 
are far apart. However, once the two objects are in the vicinity of each other, spa- 
tial decomposition techniques based on subdivision are used to solve the interference 
problem. This can become rather slow in practice as highlighted by Hahn in [21]. 

In computer animation, algorithms for polyhedral objects of complexity O(n^m^), 
where m is the number of polyhedra with n vertices per polyhedron, are described by 
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Moore and Wilhelms in [33]. For convex polytopes, linear time algorithms based on 
linear programming by Megiddo and Seidel [32, 39] and tracking closest points between 
incremental motion by Gilbert and etc. [20] have been proposed. However, they do not 
make use of coherence in their methods between successive instances. More recently, 
geometric coherence has been utilized to devise algorithms based on local features 
for convex polytopes by Baraff and Lin [3, 27]. This hcis significantly improved the 
performance of interference detection algorithms in dynamic environments. As far 
as curved models are concerned, algorithms based on interval arithmetic have been 
proposed by Duff, Herzen, Snyder, etc. [15, 22, 17]. They are relatively slow in 
practice and expect that the motion can be expressed as a closed form function of 
time. However, no practical and efficient interference detection algorithms are known 
for general geometric models in dynamic environments. 

In this paper, we present efficient algorithms for contact determination and inter- 
ference detection between two or more objects undergoing rigid motion. The models 
include polyhedral models and curved objects whose boundary can be represented 
in terms of algebraic sets. This includes algebraic surfaces and NURBS models fre- 
quently used in computer graphics and geometric modeling. No assumption is made 
on the motion of the objects. At each instance, their position is described using a 
transformation matrix with respect to the origin. For convex polytopes, the algorithm 
keeps track of closest features based on convexity of polyhedra as described in [27] 
and [28]. Concave poly topes are decomposed into convex poly topes and represented 
by sub-part hierarchical tree. 

For curved models we use a combination of hierarchical representations, algorithms 
for polytopes, and global and local methods for solving systems of polynomial equa- 
tions. In particular, we use the control polytopes of the objects and their sub-patches 
along with sub-part hierarchical tree representation to describe them for top-level 
collision detection. Once the control polytopes come into contact with each other, 
we algebraically formulate the problem of contact determination and interference de- 
tection. The problem is reduced to solving a system of overconstrained algebraic 
equations. We initially use global algebraic methods based on resultants and matrix 
computations to compute the solutions in the corresponding domain of interest and 
update these solutions using local numerical methods as the objects undergo motion. 

One of the main contribution of this paper is to present a methodology which cap- 
tures the temporal and spatial coherence between successive checks to achieve real-time 
collision detection for polyhedral objects and to determine the contact information be- 
tween general curved models efficiently using the algebraic analog of coherence. When 
the objects are not in contact with each other, the resulting algorithm for a pairwise 
test essentially runs in 6{\ -f- /(M)), where /(M) is a function of relative motion M 
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(relative homogeneous transformation matrix to describe the relative motion) which 
the objects undergo. 

We also analyze the complexity of our approaches to curved, models. The com- 
plexity of global methods is O(M^), where M corresponds to the algebraic complexity 
of the resulting system. It is a function of the degrees of the polynomials used to 
represent curved models. At each iteration the complexity of the local methods is a 
function of the algebraic degrees of the boundary surfaces and the motion undergone 
by the object pairs. 

The rest of the paper is organized as follows. We analyze the problem of interfer- 
ence detection, review literature on object models and solutions of algebraic systems 
in Section 2. In Section 3, we review the expected constant time algorithm for tracking 
closest features between convex polytopes and show how it can be extended to inter- 
ference detection between general polyhedral models. The problem of closest features 
and interference between curved models are analyzed in Section 4. In Section 5, we 
present algorithms between curved models using equation solving approaches. Finally, 
we discuss our implementation and performance in Section 6. 

2 Background 

2,1 Coherence for Detection Problem 

The simplest algorithms for interference detection are based on bounding volumes. 
Typically they correspond to bounding boxes, spheres or octrees etc. In case, two 
objects are far apart, the bounding volumes approach can be used to determine that 
there is no interference. These bounding boxes overlap only if the given pair of objects 
collide or are close to each other. As a result, any good algorithm for interference de- 
tection needs to compute the spatial relationship between a pair of objects, when they 
are close to each other and perhaps moving towards each other. Since no assumptions 
are made on the motion of the objects, it is not possible to exactly compute the time at 
which the given pair of objects would collide. Therefore, in most applications the col- 
lision detection routines are being called frequently. Due to the nature of application, 
there exists spatial and temporal coherence between successive instances. Different 
algorithms in the literature have utilized the coherence. Baraff uses cached edges and 
faces to find a separating plane between two convex polytopes [3]. However, this may 
not work when the closest features between the objects change or when there is a large 
abrupt motion. In [27], an algorithm for keeping track of closest features between con- 
vex polytopes using Voronoi regions of the features and local geometry information 
is used to determine the minimum separation between them, thus to detect possible 
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collision. This works very well for convex polytopes. However, it is difficult to apply 
this approach directly to non-convex polygonized curved models as there may be more 
than pair of closest features and the closest pair features may correspond to curves 
on the surface. In Section 5, we make use of the coherence along with local numerical 
methods and global techniques for solving algebraic systems to check whether there is 
a contact or interference between such curved models. 



2.2 Object Models 

Polyhedral models are represented using boundary representation [23]. A polyhedron 
consists of vertices, edges and faces represented in terms of parameters needed to de- 
scribe them. In addition, we also have the topological information such as adjacencies 
and incidences of all features (vertices, edges, faces). Concave poly topes are decom- 
posed into convex polytopes and represented by hierarchical sub-part trees. Curved 
models consist of NURBS surfaces and piecewise algebraic surfaces. NURBS models 
are represented in terms of control polytopes, knot vectors and order continuity [18]. 
The convex hull of the control polytopes serves as a bounding volume approximation 
for these surfaces. They are further decomposed into a series of Bezier surfaces using 
knot insertion algorithms. Each Bezier surface corresponds to a rational parametric 
surface, which can be represented in homogeneous coordinates as: 

F(»,o = (A'(»,i),K(«,0.2(«.0.>r(s.O)- (1) 

The Bezier surface is defined over the domain {s,t) € [0, 1] x [0, 1]. Algebraic surfaces 
are represented implicitly as zero sets of the form: 

f(x,y,z) = 0. 

Many commonly used models like the quadrics or torus are described using such 
formulations. 



2.3 Solutions of Algebraic Systems 

It is shown in Section 5 that the problem of interference detection corresponds to 
finding common solutions of a system of algebraic equations in a particular domain. 
The system may be a zero dimensional system with a finite number of solutions or 
an overconstrained system with no common solutions. There are different well-known 
methods for solving system of polynomials equations. Most of the algorithms in liter- 
ature are for a system of n equations in n unknowns represented as: 

Fi(a:i,X2, ...,Xn) = 0 
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( 2 ) 



F2(Xi,X2,...,X„) = 

Fn{Xi,X2,...,Xn) = 0 

Let their degrees be c?i, d. 2 , . . d„, respectively. 

Purely symbolic approaches like Grdbner bases and resultants eliminate variables 
and reduce the problem to finding roots of a univariate polynomial. These can be slow 
even for low degree systems and need exact arithmetic. On the other hand iterative 
methods like the Gauss-Newton’s method are good for local analysis and need good 
initial guesses to all the solutions. In the context of floating point arithmetic, thfe two 
main approaches are those of homotopy methods [19] and matrix computations [30]. 
The latter algorithms make use of matrix formulation of resultants of polynomial equa- 
tions and reduce the problem to computing eigenvalues and eigenvectors of generalized 
companion matrices. They have been applied to a number of geometric applications 
and perform very well in practice [30]. The complexity of the algorithm is 
where M corresponds to the algebraic complexity of the system. For dense polyno- 
mial systems, it is the Bezout bound corresponding to the product of the degrees of 
the equation; and for sparse polynomial systems, it is the BKK bound corresponding 
to the mixed volume of the Newton polytope corresponding to each equations [6, 41]. 
Homotopy methods use the solutions of a known system along with tracing algorithms 
to find the solutions of the given system. The tracing steps corresponds to Newton’s 
iteration. 



3 Collision Detection for Polyhedra 

In this section, we summarize a simple and efficient collision detection algorithm for 
convex polyhedra by tracking the closest points between them [27]. We also extend it 
to handle penetration between objects. The method works by finding and maintaining 
a pair of closest features (vertex, edge, or face) on the two convex polyhedra, in order 
to calculate the Euclidean distance between them to detect possible collision. The 
method is applicable in static environment, but is especially well suited to dynamic 
domains as the objects move in a sequence of small, discrete steps. We take advantage 
of the empirical fact that the closest features change only infrequently as the objects 
move along finely discretized paths. By preprocessing the polyhedra, the algorithm 
runs in expected constant time if the objects are not moving very swiftly. Even when 
the closest feature pair is changing rapidly, the algorithm takes only slightly longer 
(runtime is proportional to the number of feature pairs traversed which is a fuction of 
the relative motion objects undergo). 
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3.1 Preliminaries 



We represent each convex polyhedron using modified boundary representation. Each 
polyhedron data structure has fields for its features (faces, edges, and vertices) and 
Voronoi regions (defined below). 

Each feature (a vertex, an edge, or a face) is described by its geometric param- 
eters and its neighboring features, i.e. the topological information of incidences and 
adjacencies. In addition, we also precompute the Voronoi region (defined below) for 
each feature as well. 

Definition: A Voronoi region associated with each feature is a set of points closer 

to that feature than to any others. [38] 

The Voronoi regions form a partition of space outside the polyhedron according 
to the closest feature. The collection of Voronoi regions of each polyhedron is the 
generalized Voronoi diagram of the polyhedron. Note that the generalized Voronoi 
diagram of a convex polyhedron has linear size and consists of polyhedral regions. 
A cell is the data structure for a Voronoi region. It has a set of constraint planes 
which bound the Voronoi region with pointers to the neighboring cells (which share 
a constraint plane with it) in its data structure. If a point lies on a constraint plane, 
then it is equi-distant from the two features which share this constraint plane in their 
Voronoi regions. Using the geometric properties of convex sets, applicability criteria 
(explained in Sec.3.3) are established based upon the Voronoi regions. If a point P 
on object A lies inside the Voronoi region of the feature /b on object B, then /b is a 
closest feature to the point P. 

3.2 Overview 

The distance computation algorithm which tracks the closest feature pair is more fully 
described in our earlier work [27] or [28]. Here we only give a general overview of the 
algorithm. 

Our method is straightforward in its conception. We start with a candidate pair of 
features, one from each polyhedron, and check whether the closest points lie on these 
features. Since the objects (thereby the faces as well) are convex, this is a local test 
involving only the neighboring features of the current candidate features. If either 
feature fails the test, we step to a neighboring feature of one or both candidates, and 
try again. With some simple preprocessing, we can guarantee that every feature has 
a constant number of neighboring features. This is how we can verify or update the 
closest feature pair in constant time. 

When a pair of features fails the test, the new pair we choose is guaranteed to be 
closer than the old one. So when the objects move and one of the closest features 
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changes, we usually find it after one or two iterations. Even if the closest features are 
changing rapidly, say once per step along the path, our algorithm takes only slightly 
longer. In this case, the running time is proportional to the number of feature pairs 
traversed in this process. It is not more than the product of the numbers of features of 
the two polyhedra, because the Euclidean distance between feature pairs must always 
decrease when a switch is made [27], which makes cycling impossible. 

3.3 Applicability Test 

Given the overview of the algorithm, we now present a concise description of the main 
component of the algorithm — the applicability test — which is established based 
upon Voronoi region. 

A Voronoi region, associated with each feature, is bounded by the constraint planes 
derived from the feature’s neighbors. The applicability test is a simple checking process 
that verifies whether a point lies inside the Voronoi region of a given feature. With 
the preprocessing procedures to guarantee that every feature has a constant size of 
neighboring features, each applicability test is only a local test which runs in 0(1) time. 
When a point fails the applicability test of a given feature, the pointer associated with 
each constraint plane provides a new, closer feature which shares the same constraint 
plane with the previous cell. 

For a given pair of features, feat a and featg, on objects A and B, we first find 
a pair of nearest points. Pa and Fb, between these two features. Then, we need to 
verify that feats is truly the closest feature on B to Pa (i.e. Pa lies inside the Voronoi 
region of feats) and feat a is truly the closest feature on A to Ps (i.e. Ps satisfies 
the applicability test of feat a)- If either check fails, a new and closer feature (which 
is usually one of neighboring features of the previous features) is substituted, and 
the new pair of features is checked. Eventually, we must terminate with the closest 
pair, since the distance between the feature pair is strictly decreasing through each 
iteration. 

This verifying process is demonstrated in Fig.l where the previous closest feature 
candidate pair is {Fa, H)? Pa and Vi, are the two nearest points on these two given 
features. Though Pa satisfies the applicability of Vj, (i.e. Pa lies inside of H’s Voronoi 
region), Vb fails Fa’s applicability test imposed by the constraint plane CP. Therefore, 
a new candidate pair {Ea,Vb), which is closer in distance, is returned. Then, the 
algorithm is called again upon {Ea,Vb) to verify whether this is the closest feature- 
pair iteratively, until the nearest points on two features both satisfy the applicability 
tests of each other. Then, the algorithm stops and returns the closest feature-pair. 
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Object B 




Figure 1: Applicability Test: R\ and 7?2 are the Voronoi regions of Fa and Ea re- 
spectively. {Fa,Vb) —* {Ea,Vb) since H fails the applicability test imposed by the 
constraint plane CP. 
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Figure 2: Preprocessing of a vertex’s conical applicability region and a face’s cylindrical 
applicability region 

3.4 Subdivision Procedure 

In general, each face of a typical convex polyhedron has four or five edges in its 
boundary and each vertex has three or four incident edges , When a face has more 
than 5 edges in its boundary or a vertex has more than 5 incident edges, its Voronoi 
region is preprocessed by subdividing the whole region into smaller cells. In Fig. 2, Ra, 
the Voronoi region of Va which has 8 incident edges, is subdivided into i?i, R 2 , and 
R 3 ; Rb, the Voronoi region of Fb which has 8 edges in its boundary, is subdivided into 
i? 4 , i? 5 , and R^. After subdivision, each of the new cells has only 3 or 4 constraint 
planes. This subdivision procedure is a simple linear time routine which can be done 
as a precomputation step. It guarantees that when the algorithm starts, each Voronoi 
cell has a constant number of constraint planes. Consequently, each applicability 
test described above runs in constant time for updating the closest feature pairs in a 
dynamic environment. 



9 



3.5 Implementation Issues 

In order to minimize online computation time, we do all the subdivision procedures 
and build all the Voronoi regions first cts one-time computational cost. We do not need 
to reconstruct all the Voronoi regions for each polyhedron cis the objects move. We 
only need to transform the candidate features and recompute the nearest points, since 
local applicability constraints are all we need for tracking the closest pair of features. 
That is, we transform the point coordinates using relative homogeneous transformation 
matrices (available from either dynamic calculations or motion transformations as 
described in Sec 2.2) and leave constraint plane equations fixed. For example, given 
two objects moving in space, their motions with respect to the origin of the world frame 
can be characterized by the transformation matrices Ta and Tb respectively. Then, 
their relative motion can be represented by the homogeneous relative transformation 

rn rri—lrji 

This algorithm outputs the pair of closest features and also the distance between 
two objects. If the distance is less than or equal to S (a small safety margin defined 
by the user or determined by a given environment), then the objects collide. 

3.6 Hierarchical Representation for Non-Convex Objects 

We extend the collision detection algorithm for convex polyhedra to non-convex ob- 
jects using hierarchical representation. We assume that each nonconvex object is given 
as a union of convex pieces oris composed of several nonconvex subparts, each of these 
can be further represented as a union of convex subparts or a union of non-convex 
pieces. We use a sub-part hierarchy tree to represent each nonconvex object (includ- 
ing curved objects which will be discussed later Fig. 6). At each node of the tree, 
we store either a convex sub-part or the union of several convex subparts. First, we 
construct the convex hull for each node and work up the tree as part of preprocessing 
computation. We also include the convex hull of the union of sub-parts in the data 
structure. The convex hull of each node is the convex hull of the union of its children. 
For instance. The root of this sub-part hierarchy tree is the nonconvex object with its 
convex hull in its data structure. 

At each time step, we examine the possible interference by a recursive algorithm. 
The algorithm first checks for collision between two parent convex hulls. Of course, if 
there is no interference between two parents, there is no collision and the algorithm 
stops and returns the closest feature pair between two convex hulls of the objects. If 
there is a collision, then it expands their children. If there is also a collision among the 
leaves, then the algorithm recursively proceeds down the tree to determine if a collision 
actually occurs. In this recursive manner, the algorithm only signals a collision if there 
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is actually an impact among sub-parts of two objects; otherwise, there is no collision 
between the two objects. 

For complex objects, using a deep hierarchy tree with lower branching factor will 
keep down the number of nodes which need to be expanded. This approach guarantees 
that we find the earliest collision between two non-convex objects. 

3.7 Penetration Detection for Convex Polyhedra 

The core of the collision detection algorithm is built upon the concepts of Voronoi 
regions for convex polyhedra. As mentioned earlier in Sec 3.2, the Voronoi regions-form 
a partition of space outside the polyhedron. Therefore, the algorithm can possibly run 
into a cyclic loop when interpenetration occurs, if no special care is taken to prohibit 
such events. Hence, if the polyhedra can overlap, it is important that we add a module 
which detects interpenetration when it occurs as well. 

• Pseudo Voronoi Regions: 

This module can be constructed based upon similar ideas of space partitioning 
to the interior space of the convex polyhedra. The partitioning does not have to 
form the exact Voronoi regions since we are not interested in knowing the closest 
features between two interpenetrating polyhedra but only to detect such a case. A close 
approximation with simple calculation can provide the necessary tools for detecting 
overlapping. 

This is done by partitioning the interior of a polyhedron. We first calculate the 
centroid of each convex polyhedron, which is the weighted average of all the vertices, 
and then construct the constraint planes of each face to the centroid of the polyhedron. 
These interior constraint planes of a face F are the hyperplanes passing through the 
centroid and each edge in F’s boundary and the hyperplane containing the face F. 
If all the faces are equi-distant from the centroid, these hyperplanes form the exact 
Voronoi diagrams for the interior of the polyhedron. Otherwise, they will provide a 
reasonable space partition for the purpose of detecting interpenetration. 

The data structure of these interior pseudo Voronoi regions is very much like the 
exterior Voronoi regions described in Sec 3.1. Each region associated with each face 
has e-fl hyperplanes defining it, where e is the number of edges in the face’s boundary. 
Each hyperplane has a pointer directing to a neighboring region where the algorithm 
will step to next, if the constraint imposed by this hyperplane is violated. In addition, 
a type field is added in the data structure of a Voronoi cell to indicate whether it is 
an interior (pseudo) or exterior Voronoi region. 

The exact Voronoi regions for the interior space of the polyhedron can also be 
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constructed by computing all the equi-distant bisectors of all facets. But, such an 
elaborate precomputation is not necessary unless we are also interested in comput- 
ing the degree of interpenetration for constructing some type of penalty functions in 
collision response or motion planning. 

As two convex polyhedral objects move and pass through each other, the change 
of feature pairs also indicate such a motion since the pointers associated with each 
constraint plane keep track of the pseudo closest features between two objects during 
penetration phase 2 is well. However, to ensure convergence, we will need to construct 
the exact Voronoi regions of the interior space or use special case analysis to guaran- 
tee that each switch of feature necessarily decreases the distance between candidate 
features. 



4 Interference Detection for Curved Surfaces 

In this section, we analyze the problem of interference detection between curved ob- 
jects represented as spline models or piecewise algebraic surfaces. We show that these 
problems reduce to finding solutions of a system of algebraic equations. In particular, 
we present algebraic formulations corresponding to closest points determination and 
geometric contacts. 



4.1 Closest Features Formulation 



Given the homogeneous representation of two parametric surfaces, F(s,t) = (X{s,t), 
Y{s,t), Z{s,t), lV(s,t)) and G{u,v) = (X(u,u), Y{u,v), Z{u,v), iy(u,u)), the clos- 
est features correspond to points or curves on the surface. The closest features are 
characterized by the property that the corresponding surface normals are collinear. 
This can be expressed in terms of the following variables. Let 



Fii(s,t,u,w,ori) 

Fi2(-S,t,tX,U,Qri) 

F2l(5,t,U, W,Q;2) 

F22(5,t,ti,U,Qr2) 



(F(s,t) - G(tx,u)) 
ori(Gu(u,t;) X G„(u,t;)) 
(F,(s,t) X Ft{s,t)) 
Q'2(Gu(u,u) X G„(u,u)), 



where F,,F<,Gu,G„ correspond to the partial derivatives. The closest features be- 
tween the two surfaces satisfy the following equations: 



/0\ 



Fi(5,t,u,u,Qri,or2) = Fii{s,t,u,v,ai) - Fi 2 {s,t,u,v,ai) = 



0 

VV 



(3) 
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F2(5,/,u,u,ai,a2) = F2i(5,<,u,u,a2) - F22(-s,<, w,t^,0!2) = 



/ 0 \ 

0 

0 / 



This results in 6 equations in 6 unknowns. However, we are only interested in the 
solutions in the domain of interest (since each surface is defined on a subset of the 
real plane). These constraints between the closest features can also be expressed as: 



Hi(s,<,u,u) = (F(s,t) — G(u,u)) • G„(u,ti) = 0 
H2(5,f,u,t>) = (F(s,<) — G(u,u)) • G^,(u,u) = 0 
H3(s,t,u,u) = (F(s,<) -G(u,u)) •F,(s,t) = 0 

H 4 (s, f, u, t>) = (F(s, t) - G(u, u)) • F<(s, t) = 0 



( 4 ) 



where • corresponds to the dot (or inner) product. This results in four equations in 
four unknowns. 

Let us analyze the algebraic complexity of these two systems of equations cor- 
responding to closest features. Lets consider the first system corresponding to (3). 
Given 2 rational parametric surfaces F(s,t) and G(u,u), both of their numerators 
and denominators are polynomials of degree n. The numerator and denominator of 
Fii( 5, t, u, V, ai) have degree 2n and 2n due to subtraction of two rational polynomials. 
As for Fi 2 (<s, f, u, u, ai), the degrees of numerators and denominators of the partials 
are 2n — 1 and 2n respectively in the given variables (due to quotient rule). Taking 
the cross product doubles the degrees for both the numerator and denominator; there- 
fore, the degrees for the numerator and denominator of Fi2(5,f,u,u,ai) are 4n — 2 
and 4n respectively. To eliminate Oi from Fi(s,t,u,v,Qi), we get = 

cross multiplication and clearing out the de- 
nominators, we get two polynomials of degree 12n — 2 each. Once again, by the same 
reasoning as stated above for subtraction and dot product of rational polynomials, 
both the numerators and denominators of F2i(5, f , u, v, « 2 ) and F22(5, f, u, u, 02 ) have 
degrees of 4n — 2 and An. By similar method mentioned above, we can eliminate o.^ 
from F2(5, f , u, u, 0 : 2 ). We get two polynomial equations of degree 16n —4 each after 
cross multiplication. As a result the system has a Bezout bound of (12n— 2)^(16n— 4)^. 

Each equation in the second system of equations has degree 4n — 1 (obtained after 
computing the partials and subtraction of two rational polynomials) and therefore the 
overall algebraic complexity corresponding to the Bezout bound is (4n — 1)“*. Since 
the later system 4 results in a lower degree bound, in the rest of the analysis we will 
use this system. However, we are only interested in the solutions in the domain of 
interest (since each surface is defined on a subset of the real plane). 
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Figure 3: Tangential intersection and boundary intersection between two Bezier sur- 
faces 

BarafF has used the formulation in Eqn. 3 to keep track of closest points between 
closed convex surfaces [3] based on local optimization routines. The main problem is 
finding a solution to these equations (3) for the initial configuration. In general, these 
equations can have more than one solution in the associated domain (even though there 
is only one closest point pair) and the optimization routine may not converge to the 
right solution. A simple example is the formulation for the problem of interference 
detection between two spheres. There is only one pair of closest points, however 
Equations (4) have four pairs of real solutions. 

Given two algebraic surfaces, f{x,y,z) = 0 and g{x,y,z) = 0, the problem of 
closest features determination can be reduced to finding roots of the following system 
of 6 algebraic equations: 



f{xi,yi,Zi) = 0 

g{x2,y2,Z2) = 0 



^ /r(a:i,yi,2i) '' 




^ 9x{x2,y2-,Z2) '' 




= oi 


9y{x2,y2,Z2) 


< fz{xi,yuZi) y 




\ 9z{x2,y2,Z2) j 
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^ Xl ^ 




^ X2 ^ 




^ 9x{x2,y2,Z2) ^ 


yi 


— 


V2 


= 02 


9y{x2,y2,Z2) 


Zl J 




Z2 J 




{ 9z(x2,y2,Z2) ) 



Given two algebraic surfaces of degree n, we can eliminate by setting = 

_ Ajji.yi .^i) j\fter cross multiplication, we have a polynomial equation of 

P|/(l^2.y2.Z2) yx(2^2.V2.Z2) r ^ f j n 

2n — 2, since each partial has degree of n — 1 and the multiplication results in the degree 
sum of 2n — 2. Similarly, to eliminate a 2 , we set . 

and the degree of the resulting polynomial equations is n. We have six quations 
after eliminating ai and ot 2 : two of degrees (2n — 2) and four of degrees n respec- 
tively (2 from eliminating 0:2 and 2 from /(xi,yi,Zi) and y(x 2 , j/ 2 » - 2 ^ 2 ))- Therefore, 
the Bezout bound of the resulting system can be as high a,s N = (2n — 2)^72“*. In 
general, if the system of equations is sparse, we can get a tight bound with Bern- 
stein bound [5]. The Bernstein bound for Eqn. 5 is n^{n^ -1- 3)(n — 1)^. Canny 
and Emiris calculate the Bernstein bounds by using sparse mixed resultant formula- 
tion [6]. For example, the Bernstein bounds* for the case of n = 2, 3,4,5, 6, 7, 8, 9 
are 28,432,2736,11200,35100,91728,210112,435456, while the Bezout bounds are 
64, 1296, 9216,40000, 129600,345744,- • • respectively. Even for small values of n, the 
bound on the number of solutions can be large, and therefore the algebraic complexity 
of computing the closest points can be extremely high, even using the exact tight 
bound. In our applications we are only interested in the real solutions to these equa- 
tions in the corresponding domain of interest. The actual number of real solutions 
may change as the two objects undergo motion and some configurations can result in 
infinite solutions (e.g. when a closest pair corresponds to a curve on each surface, as 
shown for two cylinders in Fig. 4. ). As a result, it is fairly non-trivial to keep track 
of all the closest features between objects and updating them as the objects undergo 
motion. 



4.2 Contact Formulation 

The problem of interference detection corresponds to determining whether there is 
any contact between the two objects. In particular, it is assumed that in the be- 
ginning two objects are not overlapping. As they undergo motion, we are interested 
in knowing whether there is any precise contact between the objects. There are two 
types of contacts. They are tangential contact and boundary contact. In this section, 
we formulate both of these contact determination problems in terms of a system of 

*These figures are calculated by John Canny and loannis Emiris using their code based on the 
sparse mixed resultant formulation [6]. 
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algebraic equations. In the next section, we describe how the algorithm tests for these 
conditions as the object undergoes rigid motion. 

• Tangential Intersection : This corresponds to a tangential intersection between 
the two surfaces at a geometric contact point, as in Fig. 3(a). The contact point 
lies on the interior of each surface (as opposed to being on the boundary curve) 
and the normal vectors at that point are collinear. These constraints can be 
formulated ais: 



F(5,<) = G(u,u) 

(F,(5,f) X Ft(5,<)) • G„(u,u) = 0 ■ (6) 

(F,(5,<) X Ft(s,<)) • G„(u,u) = 0 

The first vector equation corresponds to a contact between the two surfaces and 
the last two equations represent the fact that their normals are collinear. They 
are expressed as scalar triple product of the vectors. The last vector equation 
represented in terms of cross product corresponds to three scalar equations. 
We obtain 5 equations in 4 unknowns. This is an overconstrained system and 
has a solution only when the two surfaces are touching each other tangentially. 
However, we solve the problem by computing all the solutions to the first four 
equations using global methods and substitute them into the fifth equation. If 
the given equations have a common solution, than one of the solution of the first 
four equation will satisfy the fifth equation as well. For the first three equations, 
after cross multiplication we get 3 polynomial equations of degree 2n each. The 
dot product results in the addition of degrees of the numerator polynomials (each 
of these partials has numerator polynomial of degree (2n “ 1)). Therefore, we 
get a polynomial equation of degree 6n — 3 from the fourth equation. Therefore, 
the Bezout bound of the system corresponding to the first four equations is 
bounded by = (2n)^(6n — 3), where n is the parametric degree of each surface. 
Similarly for two algebraic surfaces, the problem of tangential intersection can 
be formulated as: 



^ fx{x,y,z) \ 
fy(x,y,z) 

fz{x,y,z) j 



= a 



f{x,y,z) = 0 
g{x,y,z) = 0 
^ gr{x,y,z) \ 
9 y{x,y,z) 
9z{x,y,z) ! 



( 7 ) 



In this case, we obtain 4 equations in 3 unknowns (after eliminating a) and 
these equations correspond to an overconstrained system as well. These over- 
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constrained system is solved in a manner similar to that of parametric surfaces. 
The Bezout bound for the first three equations is M = n?{2n — 2). 

• Boundary Intersection : Such intersections lie on the boundary curve of one of 
the two surfaces. Say we are given a Bezier surface, defined over the domain, 
(s,f) € [0, 1] X [0, 1], we obtain the boundary curves by substituting s or f to be 
0 or 1. The resulting problem reduces to solving the equations: 

F(s,l) = G(«,!,) (8) 

Other possible boundary intersections can be computed in a similar manner. The 
intersection points can be easily computed using global methods. An example 
has been shown in Figure 3(b) 

Two objects have a contact if one of these sets of equations, ((6) or (8)) for para- 
metric surfaces and (7) for algebraic surfaces, have a common solution in their domain. 

In a few degenerate cases, it is possible that the system of equations (6) and (8) 
have an infinite number of solutions. One such example is shown for two cylinders in 
Fig.4. In this case the geometric contact corresponds to a curve on each surface, as 
opposed to a point. These cases can be detected using resultant methods as well [30]. 

4.3 Penetration Analysis 

In most examples, it is difficult to express the motion as a closed form function of 
time. As a result, it is almost impossible for the algorithm to compute the exact time 
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to collision as at instance we only have the information about the object’s position 
and geometry. Based on the algorithm for polyhedral models and curved models 
(presented in the next section) it is possible for us to know that the two objects are 
getting closer. As a result, the main routine of the system can check for interference 
detection at shorter time steps. However, it is possible that at a given time instance 
the two objects are not colliding and at the next instance they have penetrated into 
each other. The algebraic constraints for contact analysis, based on tangential contact 
or boundary contact, are satisfied only if the two objects have a tangential or boundary 
contact. They may no longer be satisfied, when the two objects have penetrated. In 
this section, we formulate conditions for penetration between curved models and how 
they can be used for collision detection. 

Two objects are penetrating, if there is a finite area of contact between them 
(as shown in Fig. 5). In contact analysis, we have assumed that there is a point 
contact (or a higher order contact along the boundary or on the tangent plane). The 
penetration therefore results in a one-dimensional intersection curve corresponding to 
the intersection of the two surface boundaries. In this section, we show that there is 
still a real solution to the subset of equations corresponding to (6), (8) or (7). 

Lets consider boundary intersections as in Fig. 5(b). Lets assume that there was a 
contact along the boundary curve corresponding to F(s, 1) and at the given instance 
the two surfaces have penetrated. As a result, if we compute the common solutions of 
the equations: 

F(s,l) = G(u,u) (9) 

there is a real solution in the corresponding domain, which may not correspond to 
the boundary curve of G{u,v) or there are more than one distinct solutions to these 
equations. However, any small perturbation along the boundary curve F(s, 1) would 
result in a solution to these equations in the corresponding domain and based on that 
we can detect interference. 

Now, lets consider the tangential contact between two algebraic surfaces as in 
Fig. 5(a). Similar analysis would apply to the tangential intersection of piecewise 
parametric surfaces as well. To check for tangential contact, we use the equation 
solving methods to compute the solutions of the following four equations in four un- 
knowns: 



f{x,y,z) = 0 

g{x,y,z) = 0 (10) 

/ fx(x,y,z) ^ ^ / 9 x{x,y,z) \ 

\ fy{^,y,z) ) \9y{x,y,z))' 
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Figure 5: (a) penetration resulted from tangential intersection; (b) overlap resulted 
from boundary intersection 

The solutions are substituted into the fifth equation, as shown in (6) to check for 
contact analysis. The penetration information is obtained based on following theorem. 
Lets assume that the intersection between the two surfaces is non-planar. 

Theorem 4.1 Two algebraic surfaces, f{x,y,z) = 0 and g{x,y,z) = 0 have a non- 
planar interpenetration, iff there is at least one real solution of the system of equations 
( 10 ). 

Proof: Lets take the case when the two objects have no geometric contact between 
them. As a result, there is no common real solution to the first two equations, 
f{x, j/, 2 ) = 0 and g{x, y, z) = 0. 

Lets consider the case, when the objects penetrate. Moreover, the intersection 
is non-planar. As a result, the two surfaces intersect along a one dimensional space 
curve. In a few degenerate cases, the intersection curve is higher dimensional. Lets 
denote the real curve of intersection as C. C is an algebraic space curve corresponding 
to all the common solutions of f(x,y,z) = 0 and g{x,y,z) = 0. Typically, C consists 
of one closed component in space (as shown in Fig.5). In some cases, corresponding 
to non-convex models it may consist of multiple components. Eliminating a from the 
last two equations of (10) results in 

G{x,y, z) = f,,{x, y, z)gy{x, y, z) - fy{x, y, z)g^{x, y, z) = 0. 
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The normal at a point on the intersection curve is obtained by taking the cross 
product of the two surfaces at that point. In other points, the normal function is 
defined as: 





^ fx{x,y,z) ^ 




f g^{x,y,z) 


N{x,y,z) = 




X 






fz(x,y,z) j 




{ 9z{x,y,z) ! 



It follows from this definition that G(x,j/,z) corresponds to the Z-component of the 
normal vector of the intersection curve. As a result, the real solutions of the system 
of equation (10) correspond to all those points on the intersection curve where the Z- 
component of the normal vector is zero. Since C consists of non-planar closed loops, 
there are at least two such points on each component of C. As a result, the given 
system of equations have a real solution. Q.E.D. 

The penetration criterion highlighted in Theorem 4.1 assumes that the intersection 
curve is non-planar. This is typically the case for most models (unless one of them is 
a plane or if the two models correspond to quadric surface). In fact, the non-planarity 
condition for penetration can be relaxed and it is required that the curve should 
have a point where the Z-component of the normal vector is zero. This condition 
is not satisfied, if the intersection curve lies in a plane parallel to the XK-plane. 
To circumvent this case, we take three generic linear combinations of the last three 
equations highlighted in (7) and pick two of the resulting three equations along with 
/(x,y,z) = 0 and g{x^y^z) = 0 to formulate a new system. For almost all linear 
combinations the penetration condition highlighted in Theorem 4.1 is exact. Similar 
analysis can be applied to the tangential intersection of two parametric surfaces. In 
case the two models penetrate and there is no solution for the equations corresponding 
to boundary intersection in the given domain, there is a real solution to the first five 
equations highlighted in (6) in the given domain. Again this may only fail if the 
intersection curve is planar and there is no point where the Z-component of the normal 
vector vanishes. To improve the robustness of the approach we replace the last three 
equations in (6) by their generic linear combinations. 

5 Coherence for interference detection between 
curved objects 

In most dynamic environments, the closest features between two moving objects 
change infrequently between two time frames. We have used this coherence property, 
utilizing spatial as well as temporal coherence, in designing the expected constant 
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Figure 6: Hierarchical representation of a torus composed of Bezier surfaces 

time algorithm for collision detection among convex polytopes and extended to non- 
convex polytopes as well. In this section we utilize this coherence property along with 
the algebraic formulations presented in the previous section for contact determination 
between curved models. 

Given two B-spline models, we decompose them into a series of Bezier surfaces 
using knot insertion algorithm [18]. After decomposition we use a hierarchical repre- 
sentation for the curved models. The height of the resulting tree is two. Each leaf node 
corresponds to the Bezier surface. The nodes at the first level of the tree correspond 
to the convex hull of the control polytope of each Bezier surface. The root of the tree 
represents the union of the convex hull of all these polytopes. The torus is composed 
of biquadratic Bezier surfaces, shown at the leaf nodes (Fig. 6). The convex polytopes 
at the first level are the convex hull of control polytopes of each Bezier surface and 
the root is a convex hull of the union of all control points. 



21 




1 



i 



The algorithm proceeds on the hierarchical representation, as explained in Section 
3. However, each leaf node is a Bezier surface. The fact that each surface is non-convex 
implies that the closest features of the polytopes may not be a good approximation to 
the closest points on the surface. Moreover, two such surfaces can have more than one 
closest feature at any time. Given two algebraic surfaces, we initially enclose them in 
a bounding box and test the bounding boxes for collision. At the leaf nodes, we check 
the two models for precise contact, once their control poly topes collide. 

The problem of collision detection between two Bezier surfaces is solved by finding 
all the solutions to the equations (6) and (8). If the two models are described implicitly, 
the first collision corresponds to the solution of (7). A real solution in the domain 
to those equations implies a geometric collision and a precise contact between the 
models. The algebraic method based on resultants and eigenvalues is used to find all 
the solutions to the equations (6), (8) or (7) [30]. This global root finder is used when 
the control polytopes of two Bezier surfaces collide. At that instant the two surfaces 
may or may not have a geometric contact. It is possible that all the solutions to 
these equations are complex. The set of equations in (6) represents an overconstraint 
system and may have no solution in the complex domain as well. However, we apply 
the algebraic method to the first four equations in (6) and compute all the solutions. 

The total number of solutions of a system of equations is bounded by the Bezout 
bound. The resultant method computes all these solutions. As the objects move, 
we update the coefficients of these equations based on the rigid motion. We obtain 
a new set of equations corresponding to (6) and (8), whose coefficients are slightly 
different as compared to the previous set. All the roots of the new set of equations 
are updated using Gauss- Newton’s method. The previous set of roots are used as 
initial guesses. Gauss-Newton’s method works well, if the relative motion between the 
two objects is small. This approach is similar to homotopy methods for tracing paths 
to compute solutions of a given system of polynomial equations. In case the Gauss- 
Newton’s method does not converge to a solution, we may again use the problem 
in terms of an eigenvalue problem and use the previous eigenvalues as a guess and 
use inverse power iteration. This approach has been explained in detail in [30]. The 
convergence of inverse power iteration is well understood and it converges to the 
eigenvalue closest to the guess. As a result, we are performing local computations at 
each step corresponding to a few Gauss-Newton’s iterations or inverse power iterations. 
This procedure represents an algebraic analog of the geometric coherence exploited in 
the earlier section. 

Typically, if the objects have a tangential contact and as the two objects are moving 
closer to each other, the imaginary components of some of the roots start decreasing. 
This involves the use of complex arithmetic while updating the solution with Gauss- 
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Newton’s method or inverse power iterations. Finally, a real collision occurs, when the 
imaginary component of one of the roots becomes zero. A real solution of the resulting 
system implies a collision. In case it satisfies all the equations simultaneously, there 
is a tangential contact otherwise there is penetration. The objects have a boundary 
intersection, if the equations corresponding to boundary intersection have a real solu- 
tion, which is not in the domain. As the objects undergo motion, the solution moves 
and based on the roots of the given system of equations, we can determine whether 
there is a boundary contact or penetration. 

It is possible for the objects to have multiple contacts. In such cases two distinct 
sets of solution paths converge to real solutions in the corresponding domain. As a 
result, the algorithm can easily find all the multiple point of contacts. 

The performance of the algorithm is dependent on the complexity of the global 
root finder and the number of paths we need to trace. The total number of paths 
correspond to the Bezout bound of the given system and can be fairly high for a given 
system of equations. Tracing each path at each instance can be time consuming. 
The choice of paths is made based on the solutions of the global equation solver. 
In most applications, the global root finder is used when the two objects are fairly 
close to each other. This happens when the bounding boxes corresponding to control 
polytopes overlap. As a result, we only choose those solutions as the start paths 
whose imaginary components are relatively small or whose real solutions are close to 
the corresponding domain. This is a function of the degree of the resulting surfaces. 
For higher degree models, we pick initial solutions with high imaginary components 
as well. 

Example 5.1 Lets illustrate the algorithm on two tori undergoing motion and even- 
tually colliding at a point. We assume that each torus is enclosed by a bounding box. 
Such a bounding box can either be obtained using representation of the torus as a union 
of rational Bezier surfaces (taking the convex hull of associated control points) or from 
the geometric definition. In particular^ a torus can be described as a union of points 
equi-distant from a circle in space. As a result ^ a torus can be described in terms of 
a circle C in 3 space and a radius r. The bounding box is computed by enclosing the 
circle C with the square of minimum area, say S , and its height is r along each side. 
A canonical representation of the torus, assuming a circle of radius k centered at the 
origin is 

F{x, y, z) = - r = 0. 

In Figure 7, we have shown two tori moving towards each other. The initial position 
of the first torus is described using equation 

F{x, y, z) = 64 — 20x^ + — 20y^ + 2x^y^ + J/"* + 16^^ + 2x^2^ + 2y^z^ + z”* 
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and the initial configuration of the second torus corresponds to 



G(x, y, z) = 392.0625 - 291x + 84.5x2 - 12x^ + x^ + 32.by^ - 12xy^ + 2x^y^+ 
y* — 221. bz + 84x2 — 14x^2 — 14^22; -|- Sl.bz^ — I2xz^ + 2x^z^ + 2y^z^ — 142^ + z^ 



To check for contact we add the following three equations corresponding to the first 
order differentials of F{x,y,z) and G{x,y,z): 

( — 40i + 4x^ + 4xy2 4xz^ \ / —291 + 169x — 36x2 _|_ ^j.3 _ J 2 y 2 ^ 4xj/2 -i- 842 — 28. xz — 12z2 + 4xz2 

— 40y + 4x2y + 4y^ + 4j/z2 J=al 65t/ — 24xy + 4x2y + 4j/^ — 28j/z + 4j/z2 

32z + 4x2r + 4y2z + 4z^ / \ —227.5 + 84x — 14x2 _ ^^^2 ^ _ 24xz 4x2z + iifiz — 42z2 + 4z® 



To solve for the equations, we take linear combinations of these three equations, 
eliminate a and reduce the problem of contact determination to a system of 4 equations 
in 3 unknowns. As the objects undergo motion, these system of equations is updated 
in the following manner. Let the rigid motion of the first object be represented as 



^ X ^ 




^ X \ 


y 


= Rv 


y 






V ^ / 



where R\ is an orthogonal matrix corresponding to the rotation and T\ corresponds to 
the translation. This can be expressed as 



^ X \ 




^ X 


y 


= Rl( 


y 


^ } 




V ^ / 



The corresponding position of the new object is obtained by substituting for (x,y,z) 
in F(x,y,z) and obtaining the equation F{x,y,'z) — 0 corresponding to the position 
of the torus after undergoing motion. Similarly, we compute a new representation 
G{x,y,'z) == 0 for the second torus. Notice that, it undergoes a motion represented 
by o,f^d T 2 . Given these equations, we formulate the problem of contact determi- 
nation by taking their differentials and taking the linear combinations of the resulting 
equations. The solution of the previous set of equations are used as the starting guesses. 
This process is repeated until there is a contact or penetration between the two objects 
or their bounding boxes do not overlap anymore. 
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6 Implementation and Performance 

We have implemented the algorithms described in this paper and tested their per- 
formance on many cases. The collision detection routine for convex polytopes is very 
efficient and gives us a constant time performance of 50 to 70 usee (x 10^^ sec) per pair 
of objects (independent of object complexity) on SGI Indigo2 XZ. The hierarchical 
approach works well for concave objects, given a good convex decomposition of the 
object. We are able to perform collision detection in real time using our polyhedral 
collision detection algorithm. 

Using the spatial and temporal coherence, we can keep track of the closest feature 
pairs between convex curved objects and the control polytopes of non-convex curved 
objects efficiently as well. In Fig.7(a) - 7(f) (from left to right, top to bottom) we 
demonstrate the algorithm on two tori. Each torus is represented hierarchically as 
shown in Fig. 6. In particular. Fig. 7(c) and Fig. 7(d) represent the same time instant. 
In Fig. 7(c) the convex hulls of two tori collide but in Fig. 7(d) the control polytopes 
of the subparts do not touch each other. Once two control polytopes of the subparts 
also collide Fig. 7(e), we use the global methods to find the initial solutions of (6) 
and (8) and followed by Gauss-Newton method to update the solutions at each time 
instance. It involves complex arithmetic. Although the hierarchical representation is 
in terms of biquadratic Bezier patches, we use the implicit representation of a torus 
(as that results in a system with lower Bezout bound). Each torus is represented as 
an algebraic surface of degree 4. Its equation is obtained by squaring the expressions 
in 

F{x, y, z) = + y2 _ - r = 0. 

Since the torus is a closed object, there are no boundary intersections and we keep track 
of tangential intersections for geometric contacts. The resulting system of equations 
corresponding to (6) are of degrees 4,4 and 6. The Bezout bound of this system is 
96. Using resultants, the problem is reduced to finding the eigen-decomposition of a 
96 X 96 matrix [30]. The eigen-decomposition takes slightly more than a second on 
the IBM RS/6000 and all the solutions at that instant are complex. As the tori move 
towards each other, we eventually track 8 of these solutions, whose imaginary parts 
are decreasing. This uses Gauss-Newton’s method at each time instance. Finally a 
real solution is obtained corresponding to the instance shown in Fig. 7(f). 

7 Conclusion 

In this paper we have described efficient algorithms for contact determination in dy- 
namic environments. The algorithms are based on spatial and temporal coherence 
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between successive instances and applied to polyhedral models and curved objects 
whose boundaries can be described as piecewise algebraic sets. This include NURBS 
surfaces commonly used in computer graphics and geometric modeling. The resulting 
algorithm make use of hierarchical representations, along with local numerical and 
global algebraic methods for equation solving. These algorithms have been imple- 
mented and work well in practice. 

Many issues related to contact determination for general dynamic environments 
still remain open. In an environment with N objects, we want to avoid the O(N^) 
pairwise tests at each instance. Some techniques to speed up the run time based 
on uniform spatial subdivision are proposed in [42] and algorithms of complexity 
O(Nlog^N) have been highlighted for spheres in [25] and axis-aligned bounding boxes 
in [16]. However, they are suitable for static environments but do not take advantage 
of coherence in dynamic environments. Furthermore, there are no good algorithms for 
contact determination between deformable models. 
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Figure 7: (a) - (d) Collision detection algorithm applied to two tori 
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